Perturbative approach to the nonlinear saturation of the tearing mode for any current 
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Within the traditional frame of reduced MHD, a new rigorous perturbation expansion provides the 
equation ruling the nonlinear growth and saturation of the tearing mode for any current gradient. 
The small parameter is the magnetic island width w. For the first time, the final equation displays 
at once terms of order w ln(l/w) and w which have the same magnitude for practical purposes; two 
new 0(w) terms involve the current gradient. The technique is applicable to the case of an external 
forcing. The solution for a static forcing is computed explicitly and it exhibits three physical regimes. 
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Many out of equilibrium macroscopic media display bi- 
furcations leading to the build up of macroscopic struc- 
tures. In magnetized plasmas, magnetic islands are such 
structures. They can develop in the presence of a cur- 
rent inhomogeneity through the instability of the tearing 
mode [l( which produces magnetic reconnection. This 
mode corresponds to a global magnetic perturbation that 
is resonant in a spatial region where its wave-number is 
perpendicular to the magnetic field. The magnetic is- 
land occurs in this region and stretches along the wave- 
number direction. The nonlinear saturation of a tear- 
ing magnetic island is the simplest instance of magnetic 
self-organization in a plasma. Tokamak operation avoids 
the formation of such islands, since they degrade confine- 
ment. In the reversed-field pinch (RFP), the occurrence 
of several magnetic islands leads to magnetic chaos spoil- 
ing confinement, but the formation of a single magnetic 
island is desirable, since it should provide good magnetic 
flux surfaces and a laminar dynamo. Therefore, a correct 
description and understanding of the nonlinear tearing 
mode is both important for thermonuclear fusion and for 
advancing the theory of plasma self-organization. 

The nonlinear tearing mode is classically described 
by applying resistive reduced magnetohydrodynamics 
(RRMHD) to the model of a static plasma slab, in the 
limit of small dissipation 0, El- [3- 0- El ■ The magnetic 
island region is considered as a boundary layer whose 
nonlinear features are dealt with, while the outer region 
is adequately described by linear theory only. The in- 
ner and outer solutions are then matched asymptotically. 
Rutherford 0] showed that the Navier-Stokes equation 
of RRMHD reduces in the nonlinear regime to a mere 
Grad-Shafranov equation, and proved a linear growth of 
the island width w to follow the exponential growth of 
the linear regime. The island saturation was dealt with 
approximately in Refs. |3| , and rigorously in Ref . |J| , for 
a large enough current gradient in the island region, and 
in Refs. 0, El for a vanishing gradient. 

Within the above classical setting, this Letter brings 



important novelties both in the method and in the final 
results. A new rigorous perturbation expansion using the 
magnetic island width w as a small parameter is applied 
in the inner region where no assumption is made on the 
current gradient. This technique is in principle workable 
at any order, and is applicable to the case of an external 
forcing in a plasma with a velocity profile. The explicit 
solution is provided and discussed for the case of static 
forcing. Both uniform and non uniform electric field pro- 
files are considered. The final equation describing the 
nonlinear island growth and saturation displays together 
the known term of order w ln(l/ w) [J] as well as a known 
5, 6] and two new terms of order w. Both orders have 
the same importance for physical applications, and the 
first neglected term involves a factor w 2 , which makes it 
physically smaller. Furthermore, the asymptotic match- 
ing reveals a jump of magnetic flux in the inner domain. 

We use the 2D RRMHD equations in the (x, y) plane 
that is perpendicular to the magnetic field on the reso- 
nant surface 
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where ip is the magnetic flux function, ip is the electric 
potential and plays the role of the stream function, rj and 
v are the resistivity and viscosity which may be weakly 
x-dependent, J = — V 2 i/> is the current density and J eq 
its equilibrium value. [A, B] = d x Ad y B — d y Ad x B is 
the 2D Jacobian or Poisson bracket. Units are chosen 
such that J eq (0) = 1. Furthermore, following Refs. 
and 0, it is useful to rescale the electric potential as 
</? ~~ * Vo^i where 770 = 17(0) . This makes explicit the fact 
that the vorticity and the viscosity terms in Eq. ^ are 
proportional to rj 2 and r\v respectively. These two terms 
can then be dropped by assuming the island width to be 
larger than the visco-resistive and resistive lengths 0]. 
Equations Ill2l thus become 



hM = 
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M] = rj{J m -J) (4) 
J = -Aip (5) 

In the following, two different models have been con- 
sidered: model A, in which the resistivity is constant, 
r] = 1, and J cq is given as a power series in x, J cq = 
1 + ^2 i= ibiX l , and model B, in which rj is not con- 
stant, r\ = 1 + J2iZi a i x ^ but the product i]J cq is con- 
stant (uniform electric field). In this instance, the co- 
efficients &i's can be expressed in terms of the a^s (in 
particular, a\ = —b\ and a 2 — b\ — b 2 ) In both cases, 
the origin is chosen to be the location of a null of the 
equilibrium field around which a magnetic island devel- 
ops. We thus take the equilibrium flux function to be 
ipo ~ —x 2 /2 — bix 3 /6 — b 2 x 4 /12 in a neighborhood of the 
origin, where '~' is used throughout the Letter with the 
meaning " equals plus higher order terms" . 

Small island solutions of the system of Eqs 13151 are 
conveniently obtained with the technique of asymptotic 
matching. We consider the nonlinear saturation of a tear- 
ing mode with wavenumber k. Let 5 = IV'outi — ^o| 1//2 , 
where i/'outi is the dominant part of the outer solution 
incorporating harmonics and 1 of its Fourier expansion 
in ky. We look for two classes of asymptotic solutions, 
one valid in the outer region such that \x\ >> 5, and 
one valid in the inner region (the island region) \x\ ~ S. 
Matching in the overlapping region 5 << \x\ << 1, where 
both expansions are valid, then allows to determine the 
free parameters of the problem, and in particular the ex- 
pansion parameter 5. 

In view of the matching with the inner solution, it 
is convenient to introduce the stretched coordinate £ = 
x/S, the angle x = ky, and the scaled flux function C = 
—tp/S 2 . With these definitions, the outer solution is [4|: 

Cout ~ £, 2 /2-cosx~S\n6b 1 £,cosx 

+5^/6 - (E'C/2 + A'|C|/2 + htln cos X ] 

-6 2 ln<5(6 2 /2)C 2 cosx 

+«5 2 [fo 2 e 4 /12 - {£'&i£ 2 /4 + A\ |£|£/4 

+(fc 2 /2 + b 2 - b\)e + (b 2 1 /2)e In |C|} cos X ] (6) 

where A', the usual tearing mode stability index, is the 
jump of the logarithmic derivative of i/wti at x — 0, and 
£' is the sum of the right and left values of this derivative. 

We now proceed to the nonlinear description of the 
inner region. Since the saturated mode amplitude corre- 
sponds to a bifurcation ruled by the stability index A', 
our analysis will show that S depends on A'. First, we 
need to define a suitable ordering of the fields. Since we 
are interested in tearing modes in the small island limit, 
it is appropriate to use the so-called constant-?/' approx- 
imation |]J. Denoting by ip the perturbed flux function 
and by ip' its x-derivative, we assume that ip varies little 
in the island region, Sip'/ip « 1. Since the perturbed 
current is at most of order ip' /5 and ip itself is of order S 2 , 
one concludes that J = 1 to leading order in the island 



region. Physically, this means that a tearing mode island 
does not alter appreciably the equilibrium current. At ze- 
roth order in 5, Eq. shows that ( , the leading order 
of C, satisfies <9|£o = 1, which implies Co = £ 2 /2 — cosx 
upon matching with Moreover, Eq. 0] shows that ip 
is at most of order 1. We thus rewrite Eqs. 13151 for both 
models, in terms of £ = C ~ Ccb J = J — 1, <p, and of the 



variables (£, x) 

[Co, J] = -[C,J] (7) 

k6[(o,<p] = -k5[t,tp] + ri(J eq -l)-rij (8) 

d 2 C = ~S 2 k 2 d 2 x ((o + () + J (9) 



where the Poisson brackets are now taken with respect 
to (£, x)- A few remarks are now in order. First, we note 
that whereas Eas. l7l9l can be solved in powers of 6, log- 
arithmic contributions of the form 8 n (lnd) m eventually 
appear due to the matching requirement with the outer 
solution ©. Furthermore, inspection of the structure of 
the equations shows that the exponent in In S is bounded 
by the exponent in 5: m < n. We therefore make the 
following most general perturbation expansions 

C = J2 S n (ln8) m C nm (10) 

n , m , rn < n 

J = Yl S n (lnS) m j nm (11) 

n , m , rn < n 

p = (^o+ ]T <T l (ln<5)"V nm (12) 

n,m,m< 1 n 

Finally, we note that the structure of the lowest order 
linear operator occurring in Eqs. 17191 [Co, A], makes it 
convenient to write the equations in terms of the new pair 
of independent variables (Co,x)- The Poisson bracket is 
then changed from [A, B] into £L4, B], where [A, B] is 
now defined in terms of (Co, X ) an d £ is meant as a double- 
valued function of (Co,x), £ = ±[2(C + cosx)] 1/2 - All the 
equations from now on are understood in terms of these 
new variables. Co identifies magnetic flux surfaces to the 
lowest significant order in the perturbation expansion. In 
particular, Co — 1 corresponds to the separatrix. In the 
following, we will make systematic use, for any function 
A, of the identity J c d x A dx = where C^ identifies a 
(lowest order) flux surface. The term by term derivation 
of the various contributions to the series 1101121 up to 
0(S 2 ), is now outlined. 

Order S\nS. Using Eas. 17191 and matching with © 
immediately yields jn — and Cn = ~&iC cos%. 

Order 5. To this order Eq.0is simply d x jia — 0, which 
gives jio = Ji(Co)- Writing (JSJ to the same order then 
yields 

kd x ip = h- Ji/i (13) 
Integrating ljT3|) along Cq yields 

Ji = 2^i« C o / / C'dx (14) 
/ J Cc 
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where a^ = 1 for Co > 1 and a^ — otherwise. Inte- 
gration of Eq. gives 

Ci = £^ C ° j dx - Ji dx + a{ X )li + 0(x) (15) 

where a(x) and f3(x) w iH be determined later by the 
matching conditions. 

Order £ 2 (ln<5) 2 . Equations are trivially satisfied with 
a vanishing term to this order. 

Order 8 2 In 5. Equation [3 gives the equation 



J[d x ^ix + d x j 2 i = 



(16) 



whose solution is j'21 = </{Cn + -^2i(Co) where J2i(Co) has 
yet to be determined. Using this result together with (JSJ 
and ijT3")l yields 

■W£ = kd x ((ud (o ipo ~ + b i cosx/C- (17) 
Integrating l(T7|) along Cf eventually gives 



J 2 i=6f/ ^d X C'dx (18) 

Jc Co 4 / ./c <0 

Finally, we solve for 9^21 thanks to Eq. [5] 

d&i = -6iJi(Co)cos X + /° ^rfx + 7 (x) (19) 

./ — cos y S 



Order S 2 . Eq. [3 yields 



a x j 2 = J[d x Ci 



(20) 



whose solution is j 2 = ^iCi + ^2 (Co)- Proceeding exactly 
as with J21, we get the following expressions for J2 and 

_Ic Co {jM C f}+ 6i{AJi- a Co Ci} + (62 - A6 2 )c)d x 



h 1- 1 * 



(21) 



^2 = 



J1C1 



Co 



Ji^{j}^a;-fc 2 Ccosx+e(x) 

(22) 



where A = for model A, and A — 1 for model B. 

Matching. From 10 and i|15|) . the matching condition 
gives /?(x) = and 



"(x) = &i - cosx 



£' + &! In 2 



Ji cosy ' 



cos ?/ 



dj/)(23) 



This condition still leaves the A' term unmatched. In- 
deed, it is impossible to include it in Eq. 1231 otherwise 
Ci would have a jump at £- = whereas inner solutions 



must be sufficiently smooth. We conclude that the A' 
term must be matched by higher order terms in the inner 
perturbation expansion. It is then convenient to expand 
A' as A' ~ JlndAii + tfAi. 

Wc now proceed to the next order. Upon matching 
(tm with Eq. one readily obtains 7(x) = and 



A' n = - 



J21 

lei 



dx cos x dx = (24) 



where t can be computed numerically and is approxi- 
mately 1.64, which is the result already obtained in jj]. 

The next step of our calculation is the matching of 
(|22p. One proceeds as for the previous order, but the 
calculation is more lengthy. is split into a contribu- 
tion diverging for large £'s and a finite term which after 
matching with the outer solution yields 



A[ — 2 lim \ — [ dx cos x [ 

Co-* 00 I 7T J_ n J -cosx 




W\ ]dr 



(25) 



We now reintroduce time dependance and use 124(1 and 
(12*51 to give our final result in terms of the island width 
w = AS: 



-in 2A' b\ , 
rj d t w ~ — — + —wmw- 



61S' 



w (26) 



where numerical integration gives fj, w 2.2 for model A 
and \i w 2.4 for model B. The first term on the right hand 
side was already derived in Q , the second one in , and 
the b 2 w term in 0.0]. If b\ — 0, this formula predicts a 
saturation only for 62 < 0. In the opposite case, there is 
no saturation with a small island. Furthermore, we notice 
that the £' parameter enters this equation, which means 
that the contribution of order 8 may be important even 
for a small 62 if £' is large enough. Furthermore b\ £' 
may be positive, and thus destabilizing. 

EquationlJUwas obtained by matching the first Fourier 
harmonic in x- Matching the zeroth order harmonic 
at the same orders in 6 brings in interesting physics 
too. At order 6, the matching to Eq. U3 is possible 
but for a residual contribution of the form =pS7, where 
Q = j^°(Ji — bis/2x) dx — 2&i\/2/3. This is no problem, 
since one still has the freedom to add constant contri- 
butions to the flux function @ that can be different in 
each side of the outer region. Physically, this means that 
the development of a magnetic island produces a change 
of total magnetic flux of magnitude 25 3 il. At order 5 2 , 
the zeroth Fourier harmonic of c^C2 is ^fXbiQ,. In or- 
der for this term to be matched, wc must include a new 
contribution to © of the form =FAJ 4 &if2£, which is al- 
lowed since it complies with Eq. ©-(JS)). Note that this 
contribution vanishes for model A, as a consequence of 
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the additional physical constraint of total current con- 
servation. In cylindrical geometry, the jump of magnetic 
flux corresponds to a flux production in the central part 
of the plasma by the nonlinear tearing mode. This ef- 
fect was observed in a RFP 0, and contributes to the 
toroidal field reversal, because of the global toroidal flux 
conservation 

In the presence of a static forcing, one has to add a 
solution of the inhomogcncous problem, which is conve- 
niently chosen in a unique way by setting it to zero at 
x = [{J ■ Assuming that forcing is applied at the large-x 
boundary, this solution has the form ipif(x,y) — for 
x < and ipif(x,y) = Aff(x) cos ky for x > 0, where 
fix) is the solution of the linear ideal MHD equation 
that satisfies /(0) = 0, /'(0) = 1, and Af is a parameter 
uniquely determined by the forced boundary conditions 
which can be taken positive in full generality. In the 
matching region, the complete solution is the sum of the 
unforced solution JJJJ and of the forced solution ip / , where 



At 
5 2 



+ 



S 2 h 



e + \m 



cosx 



(27) 



Moreover, a phase difference <j> between the forced and 
the unforced solutions can be taken into account by sub- 
stituting cos x with cos (x + <t>) in © ■ Matching is then 
done separately on the cosx and sin% components. The 
outcome is a set of coupled evolution equations for the is- 
land width and for the phase. One finds that fixed points 
occur for <f> — and <j) = it. For = 0, there is only one 
fixed point which is stable and whose island width is given 
by the steady state solution of l(2rj|) with the substitution 
A' -> A' + 16 A f /w 2 and £' -> £' + 16A f /w 2 . Inspec- 
tion of this equation, treating \ogw as a constant, shows 
that one can distinguish three different regimes according 

1/3 

to the value of A'. If A' A J , the saturated island 
width is essentially the one given by the unforced case, 
i.e. We oc A', with a modest increment due to the forc- 



ing. Near tearing mode marginal stability, |A'| <C A 

.1/3 



1/3 
/ ' 



w s oc Af . The resulting island is thus much larger 
than would occur without an exterior perturbation. Fi- 

1/3 

nally, in the strongly stable case, A' -C — A J , one finds 
w s cx (— Af/A') 1 / 2 , previously given in Ref. 9|. As re- 
gards the <p = 7r fixed points, the island width is given 
by the steady state solution of Ij26(l with the substitution 
A' — > A' — l6A f /w 2 and £'->£'- 16A f /w 2 . The cor- 
responding equation has physically acceptable (w > 0) 
solutions, in the number of two, only if A' exceeds a pos- 

1/3 

itive critical value of order Aj . These two fixed points 
are both unstable, the one with the largest width being 
a saddle point, unstable in the direction of the phase. 
Therefore, the island chain position always adjusts to be 
in phase with the external perturbation. 

Equation EEl has been cross-checked in two different 
ways. The analytical results have been obtained inde- 
pendently by two variants of Thyagaraja's tech- 



nique |4J . The numerical coefficients have been computed 
independently in |ll| . 

As a conclusion, we have tackled the problem of non- 
linear tearing mode saturation by using a new rigorous 
perturbation expansion. Equation 1261 brings for the first 
time the correct expression up to terms of order w 2 for 
the island width evolution. A magnetic flux jump cor- 
responding to a solenoidal effect in cylindrical geometry 
has been exhibited. The problem of the static forcing of a 
static plasma has been solved and displays three physical 
regimes. Preliminary calculations indicate that our per- 
turbation technique works also in the case of a velocity 
profile in the plasma, and of a forcing rotating at a possi- 
bly different speed than the resonant plasma region. Our 
approach can also easily be adapted to the case of cylin- 
drical geometry and this will be presented elsewhere. It is 
important to appreciate that establishing solid analytic 
techniques to solve the classic nonlinear tearing prob- 
lem opens up new routes to deal with problems whose 
modeling goes beyond conventional RRMHD. In partic- 
ular, it would be interesting to revisit neoclassical and 
two-fluid diamagnetic effects on the generalized Ruther- 
ford equation. Last but not least, our approach was not 
based on the specific J = J(ip) property but rather on a 
straightforward perturbation expansion, and may there- 
fore be applicable to more complex physical models in 
which [-0, J] ^ 0, which will be done elsewhere. 
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